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Abstract 

In this work we introduce a new approach to Dynamical Monte Carlo methods to simulate 
markovian processes. We apply this approach to formulate and study an epidemic generalized 
SIRS model. The results are in excellent agreement with the fourth order Runge-Kutta method in a 
region of deterministic solution. Introducing local stochastic interactions, the Runge-Kutta method 
is no longer applicable. Thus, we solve the system described by a set of stochastic differential 
equations by a Dynamical Monte Carlo technique and check the solutions self-consistently with 
a stochastic version of the Euler method. We also analyzed the results under the herd-immunity 
concept. 
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I. INTRODUCTION 



Epidemic systems have been systematically and mathematically formulated in a 
continuous-deterministic approach, taking immediate advantages of many numerical meth- 
ods and techniques developed to solve differential equations. The stochastic framework is 
more complex to analyze because of the required detail; therefore it could be traditionalfy 
less preferable than the deterministic ones, even being more realistic in principle 3, IJ- 
However, improved machine technology has spread the use of computationally intensive 
methods to solve a great diversity of epidemic models, and simulation techniques, as Monte 
Carlo (MC) 0, 0|, are becoming more popular in this matter. Some MC studies hide the 
effective role of the time, on time-de^ndent phenomena, reporting its results as function of 
integral Monte Carlo steps (MCS) [6]. Ambiguities of the relationship between MC time 
and real time preclude rigorous comparison of simulated results between theory and exper- 
iment. However, in the past few years, the idea of use MC methods to simulate dynamical 
processes has advanced in many publications 3, 

The aim of the present work is to present a Dynamical Monte Carlo {DMC) method for 
simulation of markovian processes. Another purpose is to incorporate explicit spatial com- 
ponents into epidemic models and analyze the dynamics of infections spread based on this 
method. We will apply the method to the compartmental Susceptible-Infected-Recovered 
(SIR) model. By the inclusion of a reflux of susceptible into the system we obtained a 
variant model: SIRS; i.e., once recovered the individual can turns back again to the class 
of susceptibles. Mean field and local interactions will be considered. We compare the re- 
sults obtained by DMC, for mean field models, with Runge-Kutta method. In cases where 
Runge-Kutta method is not applicable, the MC space-dependent results are checked self- 
consistently using a stochastic version of the Euler method 11[ and analyzed under the 
herd-immunity concept |ll2l |. 

We subdivided the work in the following way: method description (section II), Monte 
Carlo Simulation technique (section HI), epidemic models (section IV), and finally we apply 
the methodology for the solution of the SIRS model (section V). 
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II. THE METHOD 



Stochastic process approaches could simulate non-equilibrium systems, even the deter- 
ministic ones, introducing random variables to describe them in a microscopic scale. The 
macroscopic behavior of some system is resulting from averages of its microscopic proper- 
ties. Here, we will simulate systems only with markovian processes. Thus, we describe the 
evolution of the distribution of probabilities with the Master Equation: 

dm) 



dt 

J J 

where Pi{t) is the probability to find the system in the state i in the time t, and Wi^j is 
the probability of transition per unit of time. The first term on the right side describes 
the rate of all transitions to reach the considered state (increasing its probability), and the 
second term describes the rate of all transitions leaving the considered state (decreasing its 
probability). Considering Tjj as the probability of transition from the state i to j, we can 
write Wij = ^ where Xj is the characteristic time constant {lifetime) of the state i. The 
probabilities, Pi{t) and Tjj, obey the normalization conditions: Pi{t) = 1 and Tij = 1. 

We now start by choosing a convenient physical extensive microscopic quantity Ai, which 
depends only of the system's state i. Since the time must change for every successful event, 
for our purposes here, from now on we will consider only counting events related quantities. 
The mean value for a given quantity at the time t is 

A{t) = {A) = Ym)A. (2) 

i 

This equation represents the macroscopic physical quantity A{t). Differentiating both sides 

of the equation above with respect to t, we obtain 

dAjt) ^ dP^j^ .3^ 
dt ^ dt ^ ^ 

i 

and by substituting (0) in Q follows: 

^ = E E -EE (4) 

i j i j 

Defining AAij = A^ — Aj , and as i and j sweep all possible states of the system, we may 
rewrite (jlj) as 

dA{t) 



dt 



5^5^u;,^,P,AA,,. (5) 



Consider now a discretized system with interacting elements. Each element has g 
degrees of freedom given by the set {7^} with g dynamic variables 7^, i = l,2,..,g. By 
doing an element move, i.e., changing some 7^ value of a chosen element, the system will 
reach a next-neighbor microscopic state. Suppose that we are in a time scale where only 
one event occurs, and each event produces only one element move. In another words, we 
are neglecting transitions between states that need more than one element move to take 
place. Thus, let us measure "distances" among the states, say with the amount |Ay4jj|, with 
non null minimum value |AAjj| = q that defines the distance between first neighbor states. 
With the above considerations, an approach to the equation © can be done as 

dA{t) 



^Wj^iPjqdij, (6) 



dt 

in which the symbol (ij) denotes a pair of first neighbor states, and Sij = AAij/ \ AAij\. 
Now consider other physical quantity A"^ as the source for the quantity A. The use of the 
term source here is in the following sense: increasing A by the quantity q, A^ decreases by 
the same quantity and vice-versa. Thus, we can rewrite © as: 

dA{t) 



J2'tP,A]-J2'jP,A,, (7) 



dt 
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where the rate = {wj-,i). results from the average of the transition probabilities per unit 
of time, over the ensemble of first neighbor states i of j in the time t, i.e., the mesoscopic 
rates. The word ensemble here means a group of accessible configurations in a small time 
interval around the time t. We are using the time dependent ergodicity idea[l^, and in this 
sense, usually, the systems are non ergodic in non equilibrium states. The superscripts 
and "— " label the contributions to increase and to decrease the quantity A{t), respectively. 
In the particular cases where = r+ and rj = r~ are constant (or only function of the 
time) we have: 

— = r+A^ - r'A, (8) 
dt ^ ' 

which is similar to the kinetic equation for chemical reactions of first order A'^ ^ A, being 
A'^ and A the respective concentrations of the chemical elements A"^ and A. The system is in 
equilibrium when the balance at macroscopic level: r^A"^ = r~A, is satisfied. We can reach 
the equilibrium imposing the detailed balance, although this assumption is not necessary^]; 
as we will see, the equilibrium is consequence from the chosen hierarchy, that solves the 
Master Equation in any instant. 



To solve the equation (jHI) we write it in the integral form: 

A(t)-A{to)= [ Ewj^iPj{t)q6ijdt. (9) 

Jtn (ii) 



'to (ij) 

Discretizing the equation Q, we can write: 

n 

A{t) - A{to) Wj^iPjitk)qSijAtk. (10) 

fc=0 (ij) 

Let the group of possible probabilities of transition per unit of time Wj^i represented by 
the set Vt = {wj^i}, being the i and j states occurring around an instant t, with w"^^^ = 
supVt, that is, the largest probability in Vt- Each element of the time in the equation (|Tn|l 
could be 

Atk = — - — . (11) 
Starting from some initial condition, we can do the following iterative process: 

A{tk+i) = A{tk) + ^Wj^iPj{tu)q5,jAtk. (12) 

At each step k a time interval At^ is calculated using (jlip . The probabilities per unit of 
time Wj^i G Vt^. are randomly drawn using the hierarchy described in the next section of 
this work. Repeating the procedure, in a sufficient number, to get a good sample of A{tk)i 
we estimate the averages of A{tk) for every t^. This procedure is a stochastic version of the 
Euler method. 



III. MONTE CARLO SIMULATION 



In a dynamical interpretation, the MC method provides a numerical solution to the 
Master Equation. In order to do this, a sequence of events is generated based on the 
transition probabilities. The task of the MC algorithm is to create a chronological sequence 
of the distinct events separated by certain interevent times. In according to the hypothesis 
that leads to the equation (jHI), these interevent times are on a scale at which no two events 
occur simultaneously. 

To find a hierarchy to the MC algorithm, we consider n = IN, with / sweepings on the 
discretized system phase space, in which we are measuring a physical quantity represented 
in the equation ()1U|) ; in the limit of ^ oo we have the exact solution of the integral (jH)) 
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for a given initial condition. With this consideration, and using the expression the 
equation (fTUj) goes to the form: 

A{t) - Ait,) = f^Yl f Si) f ^) P^(t^)l^^r (13) 
We can, thus, create a hierarchical process choosing the transition probabilities as: 

-'j^* ...max' V-^^// 

which reproduces the correct frequencies of events in every time tk to solve ()13|) . 

To execute the MC procedure, an element is randomly selected with probability 
and thus an attempted move, with probability given by (jl4|) . is done to change the state 
from j to i. Therefore, an event that changes a dynamic variable 7^, of the chosen element, 
controls the microscopic transition j —>■ i. When the system has "degeneracy" as for the 
events occurrence, we need to decide what event will have chance (given by ()14|)) to take 



place; thus, we chose one of them with equal a priori probability [15|, supposing a local 
equilibrium over the time. The local equilibrium hypothesis means that we can measure the 
properties of the system at any instant t. Repeating the procedure, the space is swept / 
times, with the increment of time in each MCS given by (fTT|) . up to the system reach some 
desired final time. We denoted 1 MCS as a single trial to change the state of the system. 
Beginning with the same initial condition for the physical quantities, repeating the whole 
process described above, we obtain the average quantity A{t) over each instant t. As a given 
state is chosen with its correct probability in a given time, an ideal MC procedure leads to 

m / 1 \ 

A(t) - AM = Y: {{r'A^),. - (r'A),^) (^) , (15) 

where the averages are done over the ensemble of the jk states in each instant t^. 

Observing some points is necessary: first, generally different runs give different time 
results tk at the same MCS k, and we obtain the sample average with either linear inter- 
polation or extrapolation data group, in each MC realization of the system|16i]. as we will 
describe below. Second, in a complete sweep around a time tk, the value w^"-^ should be 
approximately constant in order not to change the hierarchy and consequently the result. 
Third, as the configurations do not vary drastically in few steps, the microscopic transitions 
reproduce the mesoscopic results. 



Another approach to calculate the real time consists in estimating the interevent times 
with the following rule: 

= (16) 

jk jk 

where = and A'j^ = J^-^^ or r^^^ = r~ and = Aj^^ depending, whether the result of 
the experiment increases or decreases the quantity A. The quantity /g is an e-event factor 
dependent and it obeys the relationship = 1 (normalization condition), for each time 

tk- We note that the time given by (fTB|) represents the mean waiting time for transitions 
from a given state jk to any neighbor state i; if the microscopic state stays unchanged, the 
time also does not change. We can show that this procedure leads to the same result found 
using PHI in each MCS, observing that 

^t^=y.yJ'^)(l)^tt. (17) 



Using the equation p6|) . the normalization condition for and the definition Tj^Aj^ = 
iJ^i'^jk^i we obtain the expression (fTTj). In particular, if we choose = 0, for 

every event, except e = s, we have fg=l. Under this condition, the time between s-events 
is the waiting time. Based on this, to estimate the waiting time in a coarse-grained way, we 
may define = n^/Nk] where is the number of e— events, and Nk = J2e^e is the total 
number of events, in a time interval (arbitrary) near to some time tk- 

Note that at each MCS, the minimum quantity q is either added or subtracted from the 
resulting quantity A following the prescribed hierarchy. This procedure, in according to (jl5p 
reproduces statistically the average quantity A{t). Therefore, since we have the rates, or 
the probabilities of transition per unit of time, defining the time intervals between events in 
some scale, we construct a MC algorithm to solve the Master Equation; consequently, we 
obtain the time evolution of physical quantities of the system. 

In order to define the errors on macroscopic quantities, we will do a direct approach to 
calculate average quantities. We start supposing a local equilibrium of the system over some 
instant to- With use of appropriated transition rates, we can reach any state i with proba- 
bility Pilto); so constructing several independent markov chains generates the distribution, 
which produces an ensemble of configurations over the time to- Thus, we can use directly 
the equation (j2]) to calculate A{t) at the time to- If we chose a given state i of the system 
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(18) 



A natural choice of P*(to), under equilibrium, is -Pj*(to) = -Pj(^o), obtaining 



(19) 



L* 



where L* is the number of all possible states of the system at the time to- This result 
extends readily to any time t. The states labeled by i may be considered as virtual states 
corresponding to possible data interpolation or extrapolation. The practical procedure is 
the following: at a given MC realization of the system (experiment), in the construction 
of a trajectory, labeled by i, we may get the measurements of any appropriated physical 
quantity Ai{t) obtained by either linear extrapolation or interpolation using two consecutive 
data points. After perform L Monte Carlo experiments, at some time t, the mean value of 
A is A(t) ~ X]t=i 4^- Note that if we idealize this procedure doing L — oo, we obtain the 
complete ensemble that give-us the correct mean values of physical quantities for each time 
t. Ensuring that different experiments are independent, the error for the involved quantities 
in the process for each time t could bej^: 



where A can be, for example, in this work context, the number of infected individuals. 

IV. EPIDEMIC MODELS 

The conventional treatment of epidemic systems is formulated based on a group of com- 
partments that represents each of the possible statuses, of its elements, for which we may 
assign dynamic variable values, with rates of transfer among pairs of compartments. Mathe- 
matically this subject turns into a set of differential equations. Considering a generic system 
(population, epidemic agents, etc.) and its space distributions, the temporal and space evo- 
lution characterizes any epidemic, and in each region, the density of the elements can vary 
with the time. Under this optics we considered the epidemic phenomenon as a stochastic 
process in which one random variable is the time. This focus seeks to propitiate the incor- 
poration of more details in the study of epidemic process and to allow the analysis of more 
complex models and therefore more realistic. 




(20) 
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The SIRS model considers a population with individuals divided in three classes: S 
(susceptible individuals), I (infected) and R (recovered). The evolution of the disease occurs 
according to the outline S— s>R— s>S. Based on the equation ((7j), we formalize the SIRS 
model in a quite generic way through the following group of differential equations: 

f = T,^p,R,-T,'fp>Si. (21) 



dt 

dJ 

lit 

dR 
'di 



'^ = E'fP^S,-^,f^P,I„ (22) 



J J 



^^E'-fpA-E'-rpA. (23) 

j j 

where S, I and R are the (average) number of susceptible, infected and recovered individuals, 
respectively, over each instant t. The mesoscopic rates are rj^, r^^ and r^^, for each state 
j, from S^I , I^R and R— s>S, respectively. 

In order to reproduce the deterministic model in the reference we did the following 
restrictions: 

1) the effective increase rate of those susceptible (individuals) is directly proportional to 
the number of recovered, r^'^PjRj = mR; and consequently the recovered decrease in 
the same proportion; 

2) the effective increase rate of those infected is directly proportional to the number 
of infected and a power fi of susceptible, rf^PjSj = bS'^I/N'^; and consequently the 
susceptible decrease in the same proportion; 

3) the effective removal rate of those infected is directly proportional to the number of 

infected, rY^PjIj = al; and consequently the recovered increase in the same proportion. 

Taken these restrictions to the set of differential equations fj21l - l 23|) give: 

dS „ bSf'I 

= mR- — — (24) 
dt N^" ^ ^ 

dl bS^'I 

-r = -Tr «^ 25) 

dt Nt" ^ ^ 

dR ^ „ , s 

— = al-mR, (26) 

in which /i relates to the safety-in-numbers power The conditions: dS/dt = 0, dl /dt = 0, 
dR/dt = 0, determines the steady-state solutions; the nontrivial solution occurs for finite 
values of S^, la, and R^, viz., 

Sa = {a/bf'^N, (27) 



1 + a/m 

n, . l^^N. (29) 

Depending on the removal rate a of the infectives, infection parameter 6, and renewal m, 
there exist stable solutions around the steady state that correspond to recurrent epidemics, 
or damped (fading) recurrent waves. These variant supplies oscillatory solutions that vanish 
with the time, reaching a stationary state, in which, the number of elements in each class 
stays constant. This model is a generalization of the classical SIR system readily 
recovered from ()24l — 1 26|) by setting /i = 1 , m = 0, that gives 21]: 

f - (^°) 

^ = bSI-al, (31) 

dR ^ /„„\ 
— = al. (32) 
dt ^ ' 

The SIR class of compartmental models has several deterministic and stochastic versions. 



as the SIS and the SEIR model 



20 



22,|23[. With no inclusion of spatial variables, they are 
often considered as deterministic mean field models, based on the chemical "mass action" 
principle. 

In this work, we considered epidemic processes as a result from the action of a mean 
field and the interaction among the closest individuals (local interaction). To promote the 
infection by the contact between infected individuals and susceptibles, a stochastic term is 
added to the deterministic SIRS model. Therefore, the transition probabilities per unit of 
time became 

wr-^s = (33) 
= r ^,S^-'I + A [1 - (1 - Po)1, (34) 
wi^B. = a. (35) 

The r and A parameters balance, the global (mean field) and the local (nearest neighbors) 
variables, respectively; the relation r + A = 1 is satisfied. The parameter po is the probability 
for a susceptible to become infected due to a unique infected neighbor. Therefore, (1 — 
Po)" is the probability of no infection of a susceptible if it has n infected neighbors, thus 
1 — (1 — po)" is the probability of infection of a susceptible if it has n infected neighbors. 
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The standard infection rate b, recovery rate a, exponent n and the renewal rate m are 
positive 0{1) parameters, and S{t) + I{t) + R{t) = N, with dN/dt = 0. When the renewal 
parameter is non zero (m 7^ 0), a continuous influx of susceptible rises up into the system, 
producing oscillations in the number of elements of the populational class C ={S,I, R}. 
Therefore, fading recurrent epidemics may occur before it reaches the steady (endemic) 
state. The power fi introduces a modiflcation in the original SIR model that takes in 
account nonhomogeneous mixing of susceptible and infective. 

When only the mean fleld interaction is considered, the Runge-Kutta method is enough 
to solve the SIRS model. However, the DMC method, besides to solve systems with local 
interactions, also supplies the stochastic dynamic one. 

V. APPLICATION OF THE DMC TO THE SIRS MODEL 

In this work we consider a square lattice of N = M x M sites with M = 200. The 
initial condition for the system is set up by randomly distributing Jq infectives on the lattice 
(A^ >> Jo) and the remaining sites occupied hj Sq = N — Iq susceptibles; therefore, Rq = 0. 
The simulation develops systematically by choosing one site of the lattice at a random at a 
time. Depending on its status (susceptible, infected or recovered), a trial to go to another 
status is done through a set of transition probabilities given by (jHHl - l 35|) . properly updating 
the populational class C. If the transition is successful, the system is now in a new state, 
and so we assign a time delay to this transition. We repeat the process until the system 
reaches the steady state. 

In order to construct a hierarchical process, we set the probability of transition T^^, at 
the MCS k, for a particular event a (S^I, I— s>R, or R— s>S, in our case), in accord with 
(HD), as follows: 

T^a = Wa/Wyn^^, (36) 

where &V = {ws^i, u^i^r, w^r-^s} is the transition probability per unit of time for the 
event a, and Wmax = sup P. Thus, each particular trial is gauged according to a balance of 
rates, producing a hierarchical sequence of events. Operationally, we compare T^^ against 
a random number, < TZi < 1, taken from a uniform distribution. When TZi > T^^, we 
reject the new state; otherwise accept it and calculate an incremental random time At^ 
from (Uni), with g = 1, as follows: Atf^ = -j^, or, Atf- = -gfe, or, A^^ = where 

k k k 

11 



Tl Th Th 

fi+ = ft = fs+ = = and = R. The numbers: nf^, and 

are numbers of events that increase I, decrease / and increase S, respectively, at the time t^. 
The number Afc = J2e is the total number of events, and the rates are: rf-^ = wr^s = 
rf^ = r -^S^~^I+A < [1 — (1— po)"] >jfe, = wi^r = a. The average Oj^was estimated 
doing the sum of all local interaction terms over the system configuration jk, sweeping only 
the susceptible individuals. Overall sum and search of Wmax, was in fact done only once, at 
the beginning of the simulation, after that we updated tables. We accumulated the number 
of events to obtain the factors, in two ways: first, over each 100 MCS; thus, each new 
time interval was determinated using the former calculated factor, except in the first 100 
MCS, in which we progressively calculated the factors. Second way, we let the number of 
events progressively increase, and thus, calculated the factors at each step. As the results 
agreed for both approaches, we adopted the second one because the averages converged 
faster. 

The figures, 1 and 2, show the temporal evolution of Sit), I{t) and R{t) for SIR 
and SIRS models respectively. Continuous lines represent numerical (fourth-order Runge- 
Kutta) checking solutions, and open circles correspond to the DMC simulation. Accuracies 
of the numerical solutions were checked using the steady state exact solution and the esti- 
mate of the errors were less than 0.1%. The results, shown in these figures, with respect 
to the DMC simulation correspond to an average of 20 independent trajectories, a number 
sufficient to produce soft curves and illustrate the agreement with the checking solutions. 

We introduce now the local term, with the weight A, and the variable n as an integer 
in the interval from n = up to 8, since first and second nearest infected neighbors are 
indistinguishably considered for each susceptible. From a computational point of view, the 
main consequence of introducing space-dependent variables is that the Runge-Kutta method 
is no longer applicable to the resulting model. To check the self-consistency of the approach 
we integrate numerically (^11 — I 2H() . using the Stochastic Euler method described in Section 
II; we calculated the S, I and R quantities with iterations and chose the rates randomly 
as in the MC procedure. MC solutions were checked with this method, showing excellent 



agreement 



ll|, less then 1% of difference, results not shown. The time evolution (number of 



infective) shown in Figure 3 correspond to A = 0.1, 0.5 and 0.9, and po = 0.1. Note that 
increasing A the epidemic severity reduces. Therefore, those epidemic outbreak mechanisms 
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involving only local contacts are less efficient than those whose propagation is due to some 
wider-range mechanisms. Note also, that for larger A the second peak of the curve displaces 
significantly to the right. The establishment of a protecting shield (herd immunity effect) 
may explain this effect. Depending on local contact probability, (1 — (1 —po)"'), the size 
of the removal class interferes essentially in the infection mechanism because the number 
of infectives of the neighborhood (shield effect) determines the infective character of the 
neighborhood of one susceptible. Figure 4 illustrates graphically the shield effect. 



VI. CONCLUSIONS 



In this work we examined and applied the Dynamical Monte Carlo method to the epidemic 
SIRS model. We showed that, once established the hierarchy and the relationship between 
Monte Carlo step and real time, we simulate the dynamic aspects of the system, including 
properties out of the equilibrium. Therefore, we can use the power and the generality of 
the Monte Carlo simulation to obtain the temporal evolution of deterministic or stochastic 
systems. 

We enaphasize that, here, are not required uncorrelated events as they were in the 



reference 



17l |. The results for independent runs need to be uncorrelated, so we can properly 



use the averages obtained for each time t to represent the physical quantities of the process. 
In order to do this we use a local equilibrium hypothesis, what may be at ffist glance re- 
strictive. However we may even reduce the time observation enough, for the system does 
not have time to leave some metastable states, say order of the lifetime Tf, we can so obtain 
the averaged quantities. We can obtain a good convergence to the ideal averages, in the 
practice of the simulation, by increasing the number of observations, i.e., the number of time 
experiments. 

The system studied is sufficiently general to illustrate several aspects of the real-time 
evolution determined by Dynamical Monte Carlo simulation. 

The authors gratefully acknowledge funding support from FAPESP Grant n. 00/11635-7 
and 97/03575-0. The authors would also like to thank Drs. A. Caliri and M. C. Nonato for 
many stimulating discussions and suggestions. 
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Figure Captions 



FIG. 1. SIR model. The figure shows the evolution of the number of susceptible S, 
infective / and recovered it! with time t. The numerical values for the model parameters are 
r = 0.2, b = 0.8. There is a good agreement between the MC results (open circles) and 
solutions provided by Runge-Kutta method (line) 

FIG. 2. SIRS model. The figure shows the time evolution of the S, I, R. The parameter 
values are r = 0.2, b — 0.8, m — 0.01 and /i = 2. The error between the MC results (open 
circles) and Runge-Kutta calculations are less than 0.1%. 

FIG. 3. In this figure, it is shown the effect of spatial variables for the SIRS model: the 
A and F parameters balance the local and global variables (A + F = 1); the herd immunity 
effect increases with A and it is responsible for the displacement of the curve second peak 
to the right. 

FIG. 4. The shield effect: a snapshot of the system evolution at t ~ 40, when S = 19400 
(yellow surface), R = 19200 (black) and I = 1400 (red spots) for A = 0.9. 
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